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Abstract. In this paper we consider the problem of the limits concerning the physical information that can be 
extracted from the analysis of one or more time series (light curves) typical of astrophysical objects. On the 
basis of theoretical considerations and numerical simulations, we show that with no a priori physical model there 
are not so many possibilities to obtain interpretable results. For this reason, the practice to develop more and 
more sophisticated statistical methods of time series analysis is not very productive. Only techniques of data 
analysis developed in a specific physical context can be expected to provide useful results. The field of stochastic 
dynamics appears to be an useful framework for such an approach. In particular, it is shown that modelling the 
experimental time series by means of the stochastic differential equations (SDE) represents a valuable tool of 
analysis. For example, the use of SDE permits to make the analysis of a continuous signal independent from 
the frequency sampling with which the experimental time series have been obtained. In this respect, an efficient 
approach based on the extended Kalman-filter technique is presented. Freely downloadable software is made 
available. 
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1. Introduction 

The study of the light curves of astrophysical objects has 
always been an important tool for astronomers. The rea- 
son is simple: an effective way to get insight on the struc- 
ture of a given physical system is to study its evolution 
over time. Some examples are the reconstruction of the 
structure of the binary star systems, the understanding of 
the nature of the pulsars, and the determination of the 
sizes of the central regions of the active galactic nuclei. 
However, in spite of these remarkable successes, in many 
other situations the analysis of the light curves has not 
proved to be so useful. The reason can be understood by 
taking into account that often the time evolution of a dy- 
namical system, describable in terms of a set of generic 
physical quantities {state-variables) x{t) ^, is governed by 

^ Hereafter, vector quantities will be denoted in boldface. 



a n-dimcnsional system of differential equations (state- 
equation) with the general form 

x{t)^ f[x{t),s{t),t], (1) 

where t is the time coordinate, symbol " ' " means deriva- 
tion with respect to t, and /[■] is a n-dimensional (possibly 
non-linear) function. The m-dimensional vector s{t) repre- 
sents independent processes whose time evolution does not 
depend on x{t) such as, for example, the processes that 
take place in regions external to the system of interest. 
In general, the quantities x(t) are not directly observable 
and the experimental time series {2/(^.1^^0 ^^'^ obtained 
through the measurement- equation 

yt,^h[x{tk),tk]+et,. (2) 

Here, h[-] is a Z-dimensional (possibly non-linear) function, 
{^fcjfcLo of sampling time instants, and {etj.}^Q 
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represents the measurement errors. Usually, the number / 
of available time series is smaller than n, and often 1 = 1. 
This means that the observed signals provide information 
only on a projection of the dynamics of the system of in- 
terest. Curiously, instead of developing new methodologies 
for the analysis of the observed signals in a given physical 
context, in the past much effort has been spent in an at- 
tempt to devise more and more sophisticated techniques 
for the statistical characterization of y{tk) (e.g. AR and 
ARMA modelling, maximum entropy power spectra, ...). 
The expectation of these efforts was that such a charac- 
terization could be able to provide hints on the functional 
form of /[.]. However, the results have very often been 
disappointing since the experimental time series do not 
contain all the information necessary for such a task. 

This does not mean that the classic statistical analy- 
sis of the time series is unproductive. However, it has to 
represent only a starting point, otherwise there is the risk 
that the studies on the time evolution of the astrophysical 
objects could merely consist in a collection of data and in 
the elaboration of some generic statistical measures with 
no direct physical meaning. 

On the basis of this argument. IVio et al. I l(l992l) stress 
the necessity to start directly from model as provided 
by a given theoretical model, and to use time series as a 
test for such model. More in particular, from the consid- 
eration that many astrophysical objects show evolutions 
that are unpredictable over time, they suggest to modify 
Eq. ((TJ to 

x{t)^ f[xit),uit),wit),t,9l (3) 

with u{t) and w(t) representing deterministic and ran- 
dom processes, respectively, and then to solve them via a 
numerical approach. In other words, the functional form 
of /[•] is assumed known out of a set 6 of parameters. In 
this way, once fixed the values of 0, it is possible to ob- 
tain ^^synthetic" light curves that can be compared with 
the observed ones. The reason to add the random pro- 
cess w{t), typically a continuous Gaussian white noise^, 
is a consequence of the fact that this term represents the 
interaction of the physical system of interest with its sur- 
roundings and/or the action of complex processes that 
cannot be directly included into the model (e.g. gas tur- 
bulence). In general, such processes are characterized by 
an huge number of degrees of freedom and therefore they 
can be assumed to have a stochastic nature. In practice, 
this means to study the time evolution of a given physical 
system in the context of the so called stochastic dynamics, 
i.e., through the modelling of the observed time series by 
means of stochastic differential equations (SDE). 

Although stochastic dynamics is an approach widely 
used for the study and simulation of realistic scenar- 
ios in many fields of applied science and engineering as, 

^ NB. Hereafter, with the term noise we will indicate only the 
random processes perturbing the dynamics of a given physical 
system and not the contamination due to the measurement 
errors. 



for example, fluid dynamics, structural and mechanical 
engineeri ng, avionics, material properties, financial sci- 
ences . . . (iGhanem fc Sp anos I QQltlKloeden et al. Ill997t 
ICarcia-Okialvo fc Sancho ..1999^1 . in Astronomy it is not 
so known. However, as we will show on the basis of some 
numerical experiments, even the action of quite weak noise 
sources on simple nonlinear dynamical systems can pro- 
duce deep modifications of their behaviour over time. This 
means that in real scenarios the noise component must be 
considered as intrinsic to the physics of the systems and 
not only a secondary factor. Consequently, in many situ- 
ations, stochastic dynamics could represent the only pos- 
sibility to use the experimental time series in an effective 
way. 

In fvio et al. it is suggested that the value of 

9 has to be derived from physical considerations. Here, 
we provide some tools that permits to estimate directly 
from the data. 

In Sect. El some arguments are presented that support 
the necessity of stochastic modelling in the study of the 
physical systems, and in Sec. 01 an example in the astro- 
nomical context is provided. In Sec. ^ it is shown that 
modelling the time series through discrete models is un- 
suited for most of physical systems. Hence, in Sec. an 
approach based on SDE is suggested and some tools are 
provided for using it in modelling the time series. An ex- 
ample of application is given in Sec. |S1 Finally, the con- 
clusions and some final comments are presented in SecEl 

2. Why is stochastic modelling necessary? 

In this section we will consider the time evolution of the 
simple dynamical system 

x(t) = -x^(t) +6x^{t) ~nx{t) + 6 (4) 

suffering the influence of different kinds of noise. This is a 
nonlinear model with the particularity that the associated 
potential is characterized by two regions, both of them 
with their own stable point of equilibrium (see Fig. ^ . We 
have deliberately chosen to represent an unsophisticated 
physical situation since our aim is to show that, even in 
this simplified scenario, the fact of not considering the 
noise as a fundamental component of the dynamics makes 
essentially impossible to get insights on the characteristics 
of the system under study. 

2.1. Perturbation by an additive noise process 

Figs. ^j^,h show two realizations of the stochastic differen- 
tial equation 

x{t) = ~x^{t) + 6x'^{t) - llx{t) +6 + aw{t), (5) 

where w{t) is a zero-mean, unit-variance, and continuous 
white noise process. In both the simulations the same re- 
alization of w{t) is used but two different values, respec- 
tively, 0.1 and 0.5, have been adopted for the constant a. 
Notice that in this example w{t) acts as a simple additive 
perturbing process. 
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Fig. 1. Potential function associated to the dynamical 
model Q. 

Clearly, the two signals show a different time evolu- 
tion. The reason lies in the two points of equilibrium that 
characterize the dynamics of model Q. Indeed, the only 
possibility for the system to jump from one region to the 
other is represented by the perturbation aw{t). If such 
perturbation is strong enough, then the jumps will be fre- 
quent and the bi-stable characteristic of the system will be 
revealed even by short observed signals. Conversely, if the 
perturbation is small, then it is possible that for observing 
a single jump it is necessary to wait for a long time. From 
the statistical analysis oix{t) it is possible to obtain some 
results also in this unfavo urable situation. For example, 
the Keenan test (Kccnan1 ll985(l . a test devised for veri- 
fying the nonlinearity of the time series, is able to detect 
the nonlinear nature of the data sequence shown in Fig. [St- 
at a confidence level of 95%. However, it is superfluous to 
stress that, if no jump is detected, even the most sophis- 
ticated statistical analyses will be unable to provide more 
detailed information on the functional form of model 

2.2. Perturbation by a multiplicative noise process 

The situation worsens when x{t) affects the intensity of 
the perturbation process. This is shown by Fig. ^ that 
presents a realization of the stochastic system 

x{t) = -x^{t) + 6x^{t) - llx{t) + 6 + ax{t)w{t). (6) 

Here, the realization of process w(t) is the same used in the 
previous example, and the value of a is equal to that used 
in Fig. ISb- With respect to model lO, now the standard 
deviation of the noise component at time t, being given 
by ax{t), is not constant but depends on the value of the 
signal at the same time instant. 

A comparison of Fig. with Fig. ^ indicates that, 
although the only modification regards the noise compo- 
nent, the dynamical behaviour of x{t) has suffered deep 
changes. This is illustrated in Figs. Otjb which show the 
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Fig. 2. a), b) Time series obtained from the dynamical 
model (0 with a equal, respectively, to 0.1 and 0.5; c) 
Time series obtained from the dynamical model 10 with 
(T equal to 0.5. 



0.8 




Fig. 3. Probability density functions associated a) to the 
dynamical model © and b) to the dynamical model 



probablity density functions (PDF) associated to the two 
processes. From Fig. |3Jd it is evident that the bi-stable 
structure of model Q is no more detectable. The prob- 
lem is that signal x{t), because of the term x{t)w{t), is 
no longer able to furnish direct information on the de- 
terministic part of system ©. More than in the previous 
example, this illustrates that an approach based only on 
the statistical analysis of the experimental data can be 
quite non-informative, no matter how sophisticated the 
technique applied. 

Unfortunately, it is highly probable that a situation 
like this one constitutes a typical situation for many as- 
trophysical objects. For example, it is to be expected that 
the luminosity of a given object could depend on the quan- 
tity of gas accreting onto it, while, at the same time, it 
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is conceivable that the accretion rate is influenced by the 
energy emitted by the object itself. 

3. An astronomical example 

In order to show that neglecting the stochastic compo- 
nent of a nonlinear dynamical system can be risky also 
in case of astrophysical systems, here we give an ex- 
ample based on the pair-production instability model by 
iMoskalik fc Sikoral |l986) , that has been proposed to ex- 
plain the strongly variable emission of high-energy radia- 
tion from the central regions of the active galactic nuclei 
(AGN). Since we want to maintain the readability of the 
paper also outside the context of the physics of AGN, the 
details of this model will be not given. 

According to this model, the strongly variable emission 
in the hard X-ray electromagnetic waveband of AGN's can 
be modelled via a scenario where the gas, accreting a black 
hole residing in the central regions of these objects, suffers 
a "pair production instability" (an instability due to the 
creation of electron-positron pairs) . The time evolution of 
such a system can be formalized via the following set of 
differential equations describing, respectively, 1) pair cre- 
ation and annihilation; 2) photon production, absorption 
and escape; 3) electron/positron heating and cooling; and 
4) the proton density changes: 

dn+ 
IT 

"dT 

dUe 

~dt 

dn.p 



3n., 



\ + ttR 



- 2- 



, dn+ 
~dt 



u. 



cp 



'^cbr 



Ao- — 



(7) 
(8) 
(9) 
(10) 



where 



e/cTe, Te = electron temperature, rie, n+, 



Uj, and Up are, respectively, the densities of the elec- 
trons, positrons, photons, and protons, tt — UeRffT, ctt — 
Thomson cross-section, 7i"° = pair creation rate, n!^"" — 
pair annihilation rate, Uop = rate of energy transfer from 
protons to electrons, Ucbr — electron cooling rate, Xq — 
quantity proportional to the accretion rate of the gas sur- 
rounding the central regions of AGN's , and top = electron 
energy transfer time. In the present context, the observed 
X-ray light curves can be assumed to be proportional to 
the quantity n^{t). 

An interesting point is that, for certain values of the 
parameters, this system gives rise to periodic hard X-ray 
flares. In the recent past, this model has enjoyed a certain 
fame because of this ability. However, a serious drawback 
of the above scenario is that the accretion rate of the gas is 
supposed to be strictly constant. Of course, this is a strong 
assumption and it raises some doubts on the reliability of 
the periodic behaviour under realistic conditions. 

Figs.0Ji,b show what happens to a periodic light curve, 
obtained from the numerical solution of Eqs. l(7|l- H10(l . 
when the accretion rate Ao is perturbed by a continu- 
ous, additive white noise process with standard deviation 




Fig. 4. a) Light curve n^{t) obtained by the numerical in- 
tegration of the system of equations iffjl- ljlOl) : b) The same 
light curve when Eq. Hl()(l is substituted with Eq. (|ll(l . 



cr = 0.15 Ao. Formally, this means to substitute Eq. H1Q|I 
with: 

d^A = ^X, + aw{t)]~^, (11) 



dt 



t. 



op 



that certainly represents a much more realistic assump- 
tion. It is evident that, although the perturbation is not 
so strong, the periodic behaviour of the light curve has 
almost completely disappeared. 

The indication provided by this experiment is that the 
modelling of astrophysical systems via deterministic equa- 
tions can be misleading. Indeed, some features can be fore- 
casted but, probably, are very difficult to be actually ob- 
served in experimental data. 



4. Is discrete modelling appropriate? 

In the previous sections it has been shown that neglecting 
the noise component in the dynamics of the astrophysi- 
cal systems can be risky as concerns their expected time 
evolution. Moreover, a simple statistical approach appears 
inadequate. Hence, some modelling is necessary. In this 
respect, since the experimental time series are discrete in 
nature, one could be tempted to adopt a discrete approach 
for modelling the observed signals. Actually, this is not so 
appropriate as it could appear at flrst sight. In fact, if 
the continuous nature of the signals is not taken into ac- 
count, the results provided by any method of analysis are 
in general to be expected to depend on the sampling time. 

In order to show this point it is useful to consider a 
simple stochastic model: 



X — ~6x{t) -f awit)] 



>0, 



(12) 



where 9 and a are constants. It is not difficult to see that, 
when this system is observed at a set of discrete, evenly 
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spaced time instants, its dynamics can be described in 
terms of the discrete model 



Xk+i = axk + Wk, 



(13) 



where Xk = x{tk), and {wk}k=o realization of a dis- 

crete white noise process. From the point of view of clas- 
sical time series analysis, Eq. (|13|) represents an AR(1) 
model. This fact could suggest the possibility to obtain 
information on system (|12|l by means of the classical dis- 
crete AutoRegressive modelling. Unfortunately, this is not 
true because the relationship between 9 and a ijVio et al. I 
Il992l) . 

a = exp(-6'At), (14) 

depends on the sampling time step At. The consequence 
is that, when the sampling frequency decreases, the values 
of a goes to zero. In other words, if signal x{t) is observed 
at discrete instants, it will tend to appear as a white noise 
for At oo. 

Another problem comes out from the fact that the re- 
lationship between the parameters of the continuous and 
of the discrete models can be complex and difficult to re- 
cover. For example, the second order linear system. 



x{t) + eix{t) + 0ox{t) = w{t), 



(15) 



can be shown l|Pandit fc Wu lIlQT^I) to be equivalent to a 
discrete ARMA(2,1) model, 



where. 



Xk+l - 4)Xk - 4>2Xk-l ^Wk- IpWk- 



bi = exp(/iiAi) + exp(^2At); 
62 = exp[(^i + /i2)A<]; 



V' = — ; 

V2 



■01 =(mi - ^J■2)[^J.2& 



^[(mui+fj.2)At] 



(M2e 



Ml 6' 
/ii At 



Mie 



^2 =(/i2e^^^* 
+ ^1/^2(6 



/ii At 



/t2At-j2 



:,A'2AtN2 



r-iPl 



M2j 



Mi,Ai2 



-9, ± ^ef - 400 



(16) 

(17) 
(18) 

(19) 

(20) 

(21) 
(22) 



Apart from the above mentioned dependence on At, it 
is important to note that, differently from the classical 
ARMA models, model (|16|) has only two independent pa- 
rameters: once fix and /i2 arc known, the three parameters 
01, 02 and 4' fixed. This fact implies that the classical 
ARMA(2,1) model is inadequate to represent model (|16|l 
because it assumes the independence of the characteristic 
parameters. 



The situation becomes critical in case of nonlinear sys- 
tems since, in general, it is even not possible to infer the re- 
lationship between the parameters of the continuous mod- 
els and those of the discrete ones. The reason is easy to 
understand by considering the following one-dimensional 
model 

x{t) = n[x{t),e] + a[x{t),e]w{t). (23) 

If the sampling frequency is sufficiently high, Eq. H23(l can 
be approximated by 

Xk+i ^ Xk + iJ-[xk,9]At + Wka[xk,9]{Aty/^, (24) 



i.e., with a discrete model. 



Xk+i = h{xk) + efc. 



(25) 



where h{xk) is a function that can be related back to the 
parameter vector 0, and {ek} are independent, discrete 
Gaussian random quantities. The important point is that, 
with this model, parameters 6 can be estimated through 
a classical maximum likelihood method. Things are more 
intricate if, as it usually happens in practical applications, 
At is not so small to make approximation 1)24(1 to hold. In 
fact, although often it is still possible to rewrite Eq. H23|l 
in a form, 

Xk+i = k{xk) + rjk, (26) 

in general 6 cannot be inferred from k{xk) (see also 
lTimmerll2nn(t . This has important consequences. For ex- 
ample, according to the frequency sampling, a given time 
series can display different types of nonlinearities or even 
appears as a linear process. Moreover, in general, {77^} 
are not independent, discrete Gaussian random quantities 
even when w{t) is Gaussian. 

These considerations indicate that discrete modelling 
suffers intrinsic limitations that make it unsuited to the 
description of physical systems. Gonsequently, an ap- 
proach based on continuous models is necessary. This re- 
quires the capability to estimate the parameters of a sys- 
tem of SDE from the observed time series that, in its turn, 
requires tools that permits the numerical integration of 
this kind of equations. 

5. SDE for modelling time series 

5.1. Numerical solution of SDE 

Although the theory behind SDE is complex, the numer- 
ical integration of this kind of equation does not present 
much more problems than the integ ration of the determin- 
istic differential equations (Klocdc n et al. Ill997t iHigham I 
l200l|) . Indeed, for a general SDE, 



x{t) = a[t, x{t)]dt + b[t, x{t)]w{t), 



(27) 



the simplest integration scheme, i.e. the Euler method, is 
given by 



Xt^ + Htk^tk + cTt^Wt^ 



(28) 



Here, = fi[tk,xt^], at^ = a[tk,xtj, Atk = tk+i - tk, 
and {tk}k=Q represents a set of not necessarily equispaced 
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time instants. A pitfall of this method is that its order 
of strong convergence 7 is rather small, say 0.5. For this 
reason, for the solution of the SDE's in Figs. [3 we 
have used the Milstein scheme 

- 1 , 

'fc + 7: crt.CTu 



xu+i =a;t, + fit^Atk + at^wt,^/ Atk + ^ <Tt,a't^w^^Atk, 

(29) 

that has 7 = 1. Here, at^ = a[tk,xt^], ^J-t^ = ^[tk,xtj - 
^crtk<^tk^ and the symbol " ' " denotes differentiation with 
respect to x{t). Integration schemes with higher values of 
7 are possible, however t hey are rather more c umbersome 
to implement (e.g., see iKloeden et al. lITQQTj) . The algo- 
rithms H28|l and 12911 can be easily ex tended to deal with 
systems of SDE l)Kloeden et al. Ill997(l . 

5.2. Parameter estimation in SDE 

In the past various approaches have been suggested to 
estimate the pa rameters of SDE given discrete observa- 
tions (e.g., see iBibbv fc Sorensenlll995HTimmer I I2OOC : 



iBrandt fc Santa-Clara M2002t iDurham fc GallanTTbooa , 

and reference therein). However, most of them have se- 
rious limits in the computational burden, and/or the im- 
possibility to deal with measurement errors, and/or the 
difficulty in the numerical implementation. 

Here, we present a general and flexible approach that 
is applicable to systems in the form 

x{t) = / [x{t),u{t),t, e]+(T [u{t),t, 6] w{t); (30) 

ytfc = h[x{tk),u{tk),tk,6] -f-et,,, (31) 

where x{t) G M" is the vector of the state variables. 



u{t) e 



is the vector of known input deterministic 



variables, y^^ G W is the vector of observable variables, 
e is the vector of parameters, /[•] e K", <t[-] e M"""" 
and h[-] G W are (possibly) nonlinear functions, w{t) is 
a n-dimensional, standard, continuous white noise pro- 
cess with covariance given by the identity matrix, and 
Bt^ is the vector of the measurements errors supposed to 
be zero mean, Gaussian quantities with covariance matrix 
S{uti^,tk,0). The quantitites {et^} a.ndw{t) are supposed 
to be mutually independent for all t and tk. This model 
is not so general as model Q. However, it is of interest 
in various practical applications. Moreover, as shown in 
Appendix^ through an appropriate transformation it is 
often possible to transform the more general state model 

x(t) = / [x{t),u{t), t,e]+(T [x{t), u{t), t, e] w{t), (32) 

(i.e., with (t[-] that depends also on x{t)) to the form H3U|) . 

Given a particular model structure, the maximum like- 
lihood (ML) estimation of the unknown parameters can be 
performed by finding the parameters 6 that maximize the 
likelihood function of a given sequence of measurements, 
say yt^,yt^,...,yt^,..., y^^ . By introducing the notation 

yk = [yt^,yt^_,,---,yt,,yto]^ (33) 

the likelihood function is the joint probability density 

Lie;yN)^PiyN\e), (34) 



or equivalently: 

my^)=p{yj9)l[p{y,,\yk^,,e). 



N 



(35) 



k=l 



Here the rule p{A Ci B) ~ p{A\B)p{B) has been applied 
to form a product of conditional probability densities. 
In order to obtain an exact evaluation of the likelihood 
function, the initial probability density p{yfg\6) must be 
known and all the subsequent conditional densities must 
be determined by successively solving Kolgomor oy's for - 
ward equation and applying Bayes' rule f^J azwinski Il97(]f) . 
but this approach is computationally infeasible in practice. 
Given that in Eq. H30() the term containing cr(t) does not 
depend on x{t), a more efficient alternative can be pro- 
posed. In particular, since the dynamics of model H30fl is 
driven by Gaussian, white noise processes, it is reasonable 
to assume that, under some regularity conditions, the con- 
ditional PDFs p{yj^\yk-i, 0) can be well approximated by 
Gaussians. Since the Gaussian density is completely char- 
acterized by its mean and covariance, by introducing the 
notation 

=E{y,j:y;k-i,0}, (36) 

=V{y,jyfe_i,0}, (37) 

e*. - (38) 

where E[-] and V[-] are, respectively, the the mean and 
covariance operators, the likelihood function H35(l can be 
written in the form 



N 



L{e-yr,)=p{yt^\e)\[ 



exp 



R 



tk\tk-l 



tA{Aet[R,,\,,_,]Y'' (2^)'/^' 

For a fixed 0, the quantities e^^. and Rtf,\tf._^ can be 
computed by means of an extended Kalman filter (see 
Appendix^. Further, conditioning on y^^ and taking the 
negative logarithm l{6) = — In [L{d;yk\ytg)] gives 



(39) 



N 



;(0)cx^ ln(det[fi,,|,,_J) 



(40) 



fe=i 



The ML estimate of 9 (and optionally of y^ ) can be now 
determined by solving the nonlinear optimization problem 



9 = argmin [l{9)] . 
e 



(41) 



An estimate of the uncertainty of 9 is obtained by the 
fact that the ML-estimator is asymptotically normal with 
mean 9 and covariance S given by the lower bound of the 
Cramer- Rao inequality, i.e., 

S = H \ 



Here, the Hessian matrix H — {hij} is given by 

dH{9) 



h, 



-E 



that can be estimated with 



d^l{9) \ 



(42) 



(43) 



(44) 



R. Vio et al.: Modelling of astronomical time series 



7 



6. A worked example 

To illustrate the effective potentialities of stochastic mod- 
elling in the analysis of time series, we have considered 
a sequence of X-ray observations of low mass X-ray bi- 
nary Sco X-1 (van der Klis, priv. comm.) made with 
the Proyor tional Counter Arr ay (PCA) on the Rossi-XTE 
spacecraft ( Bradt et al. Il99,'^ . The time series used in the 
experiment contains 1000 data. Sampling is regular with a 
time step of 0.015 seconds. As shown in Fig El the time se- 
ries of this object presents a power-spectrum typical of the 
quasi-periodic objects (QPO), i.e., a broad peak superim- 
posed to a steeply decreasing continuum. If one interprets 
such a fact as due to a single driving mechanism, a possible 
model for reproducing the observed signal is 

x{t) = -ax{t) + aiw{t)] (45) 
Vk = x{tk) + Acos[ujotk + <72x{tk)] + efe. (46) 

Here, a, A, cti, and a2 are unknown constants that are to 
be estimated. The frequency luq and the variance of the 
measurement errors e/j (assumed i.i.d Gaussian) are esti- 
mated through the central position of the broad peak and 
the high frequency level in the estimated power-spectrum, 
respectively. The process x{t) is assumed to be not observ- 
able. 

The idea behind this model is that the luminosity of 
the object is determined by a driving linear stochastic pro- 
cess x(t) (e.g., the accretion rate) superimposed to a pe- 
riodic process that, however, is perturbed by x{t) itself. 
Figures [3 and show that the results obtainable with 
the methodology explained in the previous section are re- 
markably good. Of course, this does not mean that such a 
simple model corresponds to a real scenario. Actually, the 
only thing that it is possible to claim is the compatibil- 
ity of model (|45|l - H4t)|) with the observations. However, we 
stress that this is the most it can be obtained from any 
technique of data analysis. 

7. Conclusions 

The time series usually available in astronomy are able 
to characterize only a subset of the system of equations 
that describe the dynamics of the physical system under 
study. For this reason, although in principle it is always 
possible to determine a statistical model able to repro- 
duce the experimental data, without any a priori physical 
model there are not so many possibilities to obtain a re- 
liable reconstruction of the physical scenario investigated. 
In general, this means that an approach to the analysis 
of time series exclusively based on the experimental data 
will provide inconclusive results, and that the practice to 
search for more and more sophisticated statistical tech- 
niques is not very productive. In many situations, the only 
possibility to get some physical insights is to carry out the 
analysis in a well defined physical context. In this respect, 
stochastic dynamics, i.e. modelling the time series with 
stochastic differential equations, appears to be a promis- 
ing tool. A benefit of working within such a framework 



Original signal 




Simulated signal 




Fig. 5. Upper panel: original (mean subtracted) time se- 
ries of SCO-Xl; Lower panel: typical time series ob- 
tainable through the fit of model H45|) - (|45|l (see text). 
Although the time series used in the analysis contains 1000 
data, here, for easiness of comparison, only the first 300 
data are shown. Time is in unit of second. 

Power-Spectrum / Original signal 




Fig. 6. Upper panel: power-spectrum |y(i^)| of the first 
5000 data of the time series of SCO-Xl; Lower panel: 
power-spectrum |y(i')| of a typical 5000 points realiza- 
tion of the process given by the fitted model H45() - (|46() . 
The frequency v is in Nyquist units. 



is that one is forced to provide a mathematical/physical 
formalization of the starting hypotheses adopted in the 
analysis of the signals (e.g. linearity, non-linearity, type 
of non-linearity, ...). In this way, there is no risk of mis- 
understandings concerning the interpretations of the final 
results. Additionally, a more direct relation between com- 
plex physical models and the limited observational mate- 
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Appendix B: Kalman Filter for estimating 
parameters in SDE 



rial is made possible permitting a more efficient interaction 

between data and theory. 

A code, implementing some efficient numerical tools 

for modelling the experimental time series with stochas- 
,. ,.rr , • f 1 J 1 J ui f of SDEs can be estimated through the maximization of 

tic differential equations, is freely downloadable from . . . . ° . 



In Sec. 15.21 it is shown that the parameters of a system 



\http://'. 



www.imm 



dtu.dk/ c tsm^ 



the likelihood function 141|) . This requires the quantities 
Etf, and Rtk,tk-i that, however, are unknown and have to 



Acknowledgements. We thank Prof. M. van der Klis for the be estimated. Here, we propose an approach based on the 

continuous- discrete extended Kalman filter. In particular, 
if at a given time instant tk the quantities 6, Xt^\t^_-i = 
E[a;tja;tfc_i], and Pt,\t,^, = E[xt^xfjxt^_^] are fixed, 
then the procedure is based on the iterated solution of the 
following sequence of equations 
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Appendix A: A multivariate transformation 

In this section, a bijective transformation 

z{t)^^[x{t),t] (A.l) 
is proposed to transform the state-equation 
x{t) = f[x{t),u{t),t, 9] + (T[x{t),u{t),t, 9]w{t), (A.2) 

to 

z{t) ^ f[z{t),u{t),t, e] + a[u{t),t, e]w{t). (A.3) 

Here, is assumed continuously differentiable with re- 
spect to t and twice continuously differentiable with re- 
spect to x{t). The covariance matrix of w{t) is assumed 
to be equal to the identity matrix. Further assumptions 
are: 

1. all the elements in cr[x{t),u{t),t,0] are strictly 
nonzero, i.e. 



(A.4) 



2. For each i there exists only one as a function of 
one and only one state-variable where v{i) is 

different for each i, i.e.. 



1. - The output prediction equations 

Vtu\t^-i = h{xt^\t^_,,ut^,tk,e) 

2. - The innovation equation 

3. - The Kalman gain equation 

4. The updating equations 

Ptk\tk = Ptk\tk-i ~ 

5. - The prediction equations 

dxt\t 



y-.Tp 1 



dt 
dt 



^ f{xt\tk,ut,t,e), te[tk,tk+i); 



(B.l) 
(B.2) 

(B.3) 
(B.4) 

(B.5) 
(B.6) 

(B.7) 



APt\t,+Pt\t,A^ + (T(T^, tG[<fe,tfe+l), 

(B.J 



aij[x{t),u{t),t,e] ^ aij[x^(i){t),u{t),t,e]; (A. 5) 

o mi r ^- r /-L\ /-L\ -L m X- J Equations (IB.7II and 1B.8II provide the quantities 

3. Ihe functions aii\x^(i)(t),u(t),t,&\ are biiective and , _ ' , ' , ' , , , 

^1 1 1 and Pf, , , \f, that can be used to start a new iteri 

a^j [x, u{t),t, 6\ are mtegrable with respect to x. 



tk+i\tk I'iiai can oe useu lo siari a new iteration of 
the sequence. Here, cr = (T{ut^,tk,0), S — S{ut^,tk,6), 

Given these assumptions, it can be shown (jNielsen et al. I 



l200l|) that the transformation 

dx 



J 



[x{t),u{t),6] 



)(t) 



(A.6) 



l,i,j = 1,2,. ..,n, fulimis Eq. \K^. 

It is interesting to notice that, after the transforma- 
tion (jA.ip , the state-equation (jA.3|l contains the same pa- 
rameters as the original state-equation (|A.2p . Moreover, 
the measurement-equation (|31|) . 



A = 



C = 



df_ 

dxt 

dh 

dxt 



(B.9) 
(B.IO) 



^-^tfc|tfc_i.^tfc,t=tfc,e 

Initial conditions for the iteration are x* 



'^t\to = and 

't\tg = J- to, which may either be pre-specified or consid- 
ered as additional parameters to estimate. 

A pitfall of the above procedure is that the matrices 
A and C have been obtained through the linearization 
can be used in its original form since the state-equation of the functions /(•) and h{-). Therefore, the approxi- 



Vt = h [xitk),u{tk),tk,6] + et^. 



x(t) is obtainable from inverse transformation x{t) 
^-'[z{t),t]. 



mated solutions obtained by solving Eqs. (|B.7|I and HB.8|I 
may bee too crude. Moreover, the assumption of Gaussian 
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conditional densities is only likely to hold for small sam- 
ple times. To alleviate these problems, a better approxi- 
mation is obtainable through a subsampling of the time 
interval [tk,tk+i), i.e., [ tk, . . . , ij, . . . , tfe+i), and the lin- 
earization of Eqs. HB.7|) . HB.8|) at each of such subsampling 
instants. This also means that the direct numerical solu- 
tion of Eqs. HB.7|I . (jB.8|) can be avoided by applying the 
analytical solutions to the corresponding linearized pre- 
diction equations 

dxt\t ^ ^ 

= f{xt\t,,Ut,,tj,e)+A{xt -Xt^) + B{ut - wt,), 

(B.ll) 



dP. 



t\tj 



dt 



= AP 



(B.12) 



where t £ [tj,tj+i), cr = cr{ut^,tj,0), S = S{ut-,tj,0), 
and 



B 



dxt 

df_ 
dut 



, .Ut ,6 



(B.13) 
(B.14) 



The solution of Eq. ljB.12|) is given by 
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T 
s 



^scra'^^lds, (B.15) 



where Tj = tj^i — tj, and $s = e^^. The solution of 
Eq. HB.11(I is more difficult to find. If A is nonsingular, it 
is given by 

Xt^^,\t,^Xt^\t,-TsA-'Ba + A-\^,-I){A-'Ba + f), 

(B.16) 

where 



a = 



and / = f{xi\f. , Ut- , tj, 6). Things are mo re complex if A 
is sing ular. More details can be found in iKristensen et al~l 
ll2004l). 
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